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Quantum energy transfer in a chain of two-level (spin) units, connected at its ends to two thermal reser- 
voirs, is analyzed in two limits: (i) In the off-resonance regime, when the characteristic subsystem excitation 



i-^ . energy gaps are larger than the reservoirs frequencies, or the baths temperatures ai^e low. (ii) In the reso- 

O , nance regime, when the chain excitation gaps match populated bath modes. In the latter case the model 



is studied using a master equation approach, showing that the dynamics is ballistic for the particular chain 



^ ■ model explored. In the former case we analytically study the system dynamics utilizing the recently devel- 

ff^ , oped Energy-Transfer Born-Oppenheimer formalism [Phys. Rev. E 83, 051 1 14 (201 1)], demonstrating that 

cn 

■^ \ energy transfers across the chain in a superexchange (bridge assisted tunneling) mechanism, with the energy 

Q ' current decreasing exponentially with distance. This behavior is insensitive to the chain details. Since at 

low temperatures the excitation spectrum of molecular systems can be truncated to resemble a spin chain 
model, we argue that the superexchange behavior obtained here should be observed in widespread systems 

X" 

\^ • satisfying the off-resonance condition. 



I. INTRODUCTION 



The scaling of the 



energy current with system size is of interest for developing applications 



in energy conversion yj], molecular electronics [|20, and reaction dynamics 



li. 



In the context of 



biological macromolecules, understanding pathways and efficiency of heat flow is important for 
controlling signal transmission and functionality in biomolecules [|4J]. In nanoscale electric devices 
significant power dissipation may lead to the system disintegration. Designing efficient routes for 
energy transfer, away from the conducting object, is essential for a stable operation [|5D. 

Resolving the size effect of heat conduction in molecular chains has been the subject of re- 
cent experiments la lZO. For example, Wang et al. [8] has recently studied the kinetics of heat 



m. 



transfer from a metal substrate to self-assembled hydrocarbon monolayers of increasing sizes 
concluding that heat travels ballistically along the chain llioll . Vibrational energy transport in 
peptide helices was measured by employing vibrational probes as local thermometers at various 
distances from a heat source [11]. For this protein system it was concluded that heat propagates in 
a diffusive-like process. 

Considering electron transfer diCro^^ a 1-dimensional (ID) conductor (bridge) connected at the 
edges to electron reservoirs, multitude of theoretical, numerical, and experimental studies have 
demonstrated that in the resonance limit, applicable for ohmic reservoirs at high temperatures, 
(quasi) ID chains conduct electrons anywhere between a ballistic to a diffusive manner, depending 
on the details of the internal interactions. In contrast, in the off-resonance regime, when the 
electronic levels of the bridge lie high, above the populated states of the reservoirs (the donor and 
acceptor states in a donor-bridge acceptor complex), a deep tunneling mechanism should take off 



In this paper we focus on the analogous problem of energy transfer in molecular chains con- 
nected at the two ends to thermal reservoirs (baths), maintained at distinct temperatures, realized 
by solids, metals, nanoparticles [|l3n or large molecular complexes |l6|]. As a simple model for the 



central object we consider a chain of several two-level units (spin 1/2 particles). The units are 
coupled through nearest-neighbor coupling terms, assumed to be weak compared to the on-site 
energies. Similarly to the electronic case, we expect that distinct transport mechanisms will dom- 
inate at different parameter regimes. For a schematic representation of the chain model and the 
relevant excitation spectra see Fig. [B 

Before proceeding, we carefully clarify our terminology: "Resonant" regime refers here to the 
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FIG. 1: Top: Chain of two-level particles or spin- 1/2 objects coupled at the edges to heat baths maintained 
at different temperatures. Bottom: In the off -resonance regime, applicable, e.g., at low temperatures, the 
chain's modes, relevant for the transport process, match unpopulated bath modes. In the opposite resonance 
regime the chain's modes overlap with populated bath modes. The coloring of the baths' spectral func- 
tion represents thermal population, where darker color reflects larger occupation. The off-resonant model 
could be also realized considering reservoirs whose spectral functions have a low cutoff, below the chain's 
characteristic excitation gaps. 



case where the modes of the chain, responsible for the heat transfer dynamics, lie in resonance 
with the occupied baths' modes. In contrast, "Off- Resonance" conduction refers to the case where 
the occupied modes of the thermal reservoirs' are low, below the typical excitation frequencies of 
the enclosed object. Overall, there is always a conservation of energy in our system, transferred 
between the two reservoirs (donor and acceptor). Thus, irrespective of the bridge energetics, we 
always consider here a "resonance energy transfer" (RET) process lll4 ll5ll. in the sense that there 
are no energy loss mechanisms, e.g., it is a non-radiative process. 

In the resonant regime numerous theoretical and computational studies have demonstrated that 
energy transfer between two reservoirs, mediated by the excitation of the interlocated object, may 
follow a ballistic J oc N^, ohmic J oc A^~^, (or somewhere in between, J oc N"-^^, a > 0) 
mechanism. This was done in the context of vibrational energy transfer lll6n and thermal transport 



in spin chain 



191], for both classical and quantum systems, using, e.g., molecular d yna mics 



simulations 11611 . the quantum master equation method 1191-2111 . the Green-Kubo formula 122 , 
and the density matrix renormalization group method [|24ll . In the off-resonant regime Simula 



m. 



tions on purely harmonic systems indicated on a tunneling dynamics of heat transfer [llOll . In the 
presence of interactions, off-resonant quantum heat transfer dynamics has been treated by adopt- 
ing e.g., the complex machinery of the Green's function approach [1251 12611 . or by using mixed 
quantum-classical simulations MM- Typically, such methods only provide numerical results, hin- 
dering direct picture of the microscopic processes involved. Responding to this challenge, we have 
recently developed a simple analytic method for describing energy transfer in nonlinear systems in 
the off-resonant regime ll28ll . This method, an extension of the Bom-Oppenheimer (BO) principle 
[I29I1 to energy transfer problems, can treat general subsystems (impurity, chains) with intrinsic an- 
harmonicities, as well as cases where the subsystem is nonlinearly coupled to the reservoirs. The 
outcome of the method is a Landauer type expression, incorporating nonlinear interactions BOll . 
allowing for the identification of different scattering processes 112 811 . In what follows we refer to 
this method as the Energy-Transfer Born-Oppenheimer (ETBO) scheme. 

In this paper we aim in deriving scaling laws for the behavior of the energy current with size 
for ID molecular systems, primarily focusing on the off-resonant limit [l3l|]. For simplicity, we 
consider the isotropic XY spin chain and its variants as a prototype for a homogeneous and linear 
molecular chain [13211]. This is a relevant physical description since at low temperatures or in the 
off-resonance limit the energy spectra of the interlocated system can be truncated, as transport 
predominantly occurs through the lowest excitation states. Considering a spin chain between two 
thermal reservoirs, we study the energy transfer behavior in two different limits: (i) We assume an 
off-resonance scenario, and obtain the energy current adopting the ETBO approach. In this case 
we demonstrate that the energy current decays exponentially with size, a footprint of the tunneling 
mechanism, (ii) In a resonant situation we utilize a standard master equation approach and show 
that the isotropic XY chain behaves as a ballistic conductor, providing a fixed current for different 
sizes. While we present our study in the context of steady state heat transfer, the results are also 
useful for interpreting energy transfer rates in donor-bride- acceptor complexes 113 311 



The tunneling behavior of the energy current resolved in the off-resonance regime [|l4l ll5. 



m 



resembles the McConnell superexchange result [|35|1 . observed in electron transport experiments 
in numerous systems, including monolayers 1360. proteins [|370, and DNA [38l]. Since the thermal 
superexchange result does not depend on the details of the chain model, we expect it to show 
up in different physical systems at low temperatures, including molecular wires, spin chains, and 
biomolecules. 

Our study here is presented in the context of thermal energy transfer. However, the analy- 
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FIG. 2: Energy spectrum of the isolated isotropic XY spin chain with iV = 2, A'^ = 4, A^ = 6, and A^ = 8 
units. Other parameters are e = 2 and n = 0.1. The different manifolds include different numbers of 
excitations on the chain, from zero up to N. 

sis and results are valid for describing general excitation energy transfer problems (vibrational, 



electronic), in bulk-molecule-bulk junctions and donor-bridge-acceptor systems 1114 , 



15L 



SJ. Ex- 



periments on cr-bond or vr-bond bridges, connected to donor and acceptor c 



iromophores reported 
on excitation transfer rates which are exponentially decreasing with size Il39l l40ll. This behavior is 
rigorously recovered here. 

The structure of this paper is as follows. In Section II we describe the spin chain model, 
serving as a prototype for studying energy transfer and thermal conduction in linear chains. In 
Section III we study the off-resonant case using the ETBO method. Perturbative analytic results 
are supported by numerical simulations. Section IV treats the resonant limit, adopting a master 
equation approach. Section V concludes. 



II. MODEL 

Consider a small subsystem, representing e.g., a molecule, placed in between two thermal 
reservoirs (e.g., solids, large complexes) maintained each at a fixed temperature T^ = /3~^ (u = 
L,R). The total Hamiltonian is given by 



H = Hs + Hl + Hr + Vl + Vi 
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FIG. 3: Energy spectrum of the isolated anisotropic Heisenberg spin chain with A^ = 4 (a) and A^ = 10 (b). 
Panel (c) zooms on a small portion of the spectrum for A^ = 10, manifesting the gap closure with increasing 
6. Other parameters are e = 2, k = 0.1 and 6 = (o), 6 = 1 (+) and 6 = 2 (dotted). 

Hs is the Hamiltonian of the subsystem and Hi^ stands for the u heat bath. V^ couples the subsys- 



tem and the i' reservoir. In particu 



of the isotropic XY spin-chain [|32h 



ar, in what follows we focus on the thermal transport properties 



N N-1 

Here crj'^'^ are the Pauli matrices for the j spin. The first term describes the onsite spectrum at 
each site, a two level system with a spacing e. The second term provides the hopping interac- 
tion between neighboring sites with an exchange strength k. We generally assumes that k < e, 
allowing for a meaningful description of the chain in terms of its subunits. The chain is coupled 
to two independent thermal reservoirs at sites 1 and N. Our derivation below does not make any 
assumption regarding the structure of these baths. For example, they may each contain a collection 
of independent harmonic oscillators 



Hu = J2^^^lj^'''r 



(3) 



jeu 



One may also consider fermionic reservoirs, a source of electronic excitations. System-bath inter- 
actions are assumed to take the following form, 



Vl = SiBl, Si = Xlc^i 
Vr = S^Bji; Sn = Aiju;^, 



(4) 



where X„ parametrizes the system-bath coupling strength, assumed to be a real number, 5*1 (Sn) 
are subsystem operators coupled to the left (right) reservoirs, Bl {Br) are L (R) bath operators, 
which need not be specified at this stage. This form assumes that the interaction with the reservoirs 
can generate (absorb) an excitation at the leftmost or rightmost sites of the chain. 



Using the Jordan-Wigner transformation 14111 . the (isolated) isotropic XY chain can be reduced 
to describe spinless free fermions. However, our analysis is still not trivial since the chain is 
coupled, potentially in a nonlinear way, to thermal reservoirs which are not necessarily modeled 
by isotropic XY chains by themselves. Thus, the total Hamiltonian cannot be transformed into a 
collection of noninteracting fermions, and the model is generally non-integrable. Furthermore, the 
ETBO analysis can be carried out for studying transport through the ID anisotropic Heisenberg 



model 1I32I1 . 



^^ = i E -^ + i E (-f-l+i + -14m) + Y-^-l+i- (5) 

The spectrum of Hs with 5 = is exemplified in Fig. [2] for several sizes. We note that for 
e ^ K subsystem energies are grouped into manifolds, each including eigenstates with a particular 
number of excitations: At the bottom of the spectrum lies a zero excitation state with an energy 
of £"0 ~ i~^^)- Above, we identify states including a single excitation on the chain, with their 
energy centered around Eq + e. The next manifold includes the two-excitation states, and so on 
and so forth. For example, for A^ = 4 there are five manifolds with zero (bottom) to four (top) 
excitations residing on the chain. Note that given the form of the system-bath interaction operator 
[Eq. dH)], the thermal baths, the excitation resources, can only translate the subsystem between 
(neighboring) manifolds, adding or absorbing a single excitation at a time. The reservoirs cannot 
directly drive transitions within each manifold. Since for k ^ e the typical gap between manifolds 
is ~ e, this energy scale is identified as the characteristic frequency of the subsystem, controlling 
the transport properties of the model. 

In Fig. [3] we show that this picture is retained when the exchange anisotropy parameter 5 is 
relatively small, for short chains. Then, one can still identify manifolds with different number of 
excitations, as the gap between bands is larger than energy differences within each band. However, 
for large 5 and for long chains the spectrum becomes more involved, and the gaps between different 
excitation states diminish. In what follows we restrict ourselves to situations where the gaps 
between manifolds are maintained, ~ e, larger than the spacings within each band (bandwidth 
~ k). Practically, we study chains of A^ < 10 units, with large onsite gaps e ^ k and small 
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anisotropy parameter < 5 < 1. 

The spin chain serves as a prototype model for exploring energy transfer through a linear molec- 
ular junction, at low temperatures. In what follows we study the steady state behavior of this model 
in two different limits: 

(i) Off-resonance case, with e ^ tUc or T^ < e. Here cUc is the reservoirs cutoff frequency. 
In this limit energy transfer takes places between low frequency modes of the two reservoirs, 
mediated by a (high frequency) subsystem. For treating the dynamics in this scenario, we adopt 



the recently developed Energy-Transfer Born-Oppenheimer method 112811 . In Sec. Ill we show 
that in this regime heat is transferred via a coherent-superexchange mechanism, with the current 
exponentially decreasing with chain size. 

(ii) Resonance regime, with e < cOc and Tj, > e. Under these conditions baths' modes in 
resonance with the subsystem frequencies are populated, responsible for the subsystem excitation 
and relaxation processes. This resonance energy transfer process can be treated within the Born- 



Markov approximation scheme [|42l . 



4311 . In Sec. IV we show that in this case the isotropic XY 



model transfers energy in a ballistic manner. 

III. OFF-RESONANCE REGIME: ENERGY-TRANSFER BORN-OPPENHEIMER SCHEME 

A. Method 

We describe the principles of the Energy-Transfer Bom-Oppenheimer method as developed in 
Ref. [281, then apply it onto the spin-chain model, to obtain the behavior of the current as a func- 



tion of size. Generally, the BO approximation ll29ll is based on the recognition of timescale sepa 



ration. In isolated molecules, the "traditional" BO approximation relays on the mass separation of 
electrons and atomic nuclei. In this context, one assumes that the electron cloud instantly adjusts 
to changes in the nuclear configuration, and that the nuclei propagate on a single potential energy 
surface associated with a single electronic quantum state, obtained by solving the Schrodinger 
equation with fixed nuclear geometries. 

This principle can be adopted for treating quantum thermal transport in (potentially strong) 



interacting systems driven to a steady state by a temperature bias i28|l . The method is applicable 



in the off-resonant regime, where the characteristic frequencies of the impurity object are high 
relative to the cutoff frequencies of the reservoirs e ^ Wc 13 IQ . This implies a timescale separa- 



tion, as the subsystem dynamics is fast, while the bath motion is slow. The ETBO approximation 
follows two consecutive steps: First, we consider the fast variable and solve the subsystem eigen- 
problem while fixing the reservoirs configuration, to acquire a set of potential energy surfaces 
which parametrically depend on the bath coordinates. In the second step we adopt the adiabatic 
approximation and assume that the baths' coordinates, the slow variables, evolve without changes 
on the subsystem state. We then solve the energy transfer problem between the reservoirs on a 
fixed potential surface. 

In what follows we denote by q subsystem coordinates and by Q^ the u bath coordinates. These 
are collections of displacements and momenta operators. The baths operators which are coupled 
to the subsystem. By, are functions of the Qy coordinates. We also collect in Q the coordinates of 
both reservoirs. Fixing the bath coordinates, we identify the 'fast' contribution to Eq. ([T]) as 

H^{q, Q) = Hsiq) + ^ K(g, Q.), (6) 

and solve the time-independent Schrodinger equation 

iiM.Q)\9n[q.Q)) = Wn{Q)\gn{q,Q)), (7) 

to acquire a set of "potential energy surfaces", Wn{Q), and states \gn{q, Q))- Note that the po- 
tentials Wn mix the left and right system-bath interaction operators. Moreover, they are not nec- 
essarily linear in Ql and Qji. These potentials are the analogs of the electronic potential energy 
surfaces obtained in molecular structure calculations, which parametrically depend on the nuclear 
coordinates. Similarly, the reservoir coordinates Q are treated as parameters in Eq. ([7]). Assuming 
that the surfaces are well separated, we presume the adiabatic ansatz and write the total density 
matrix as 

p{t) = \9n{q,Q))pliQ,t){gn{q,Q)\, (8) 

where the bath density matrix obeys the Liouville equation 

pl{Q,t) = e-'''Boip^(0)e'H-^\ (9) 

h=l, with the effective Hamiltonian 

H^o = Hl{Ql) + Hr{Qr) + W^{Ql. Qr). (10) 

Here pb{0) = pl x Pij is a factorized initial condition with p,^ = e~^/Tr,^[e~^], the equilibrium- 
canonical distribution function of the u bath. In what follows, we assume that the baths coordinates 



evolve on the ground potential surface, simply denoted by W . The effective Hamiltonian (UO)) 
including the ground potential surface W will be similarly denoted by Hbo- For brevity, we also 
omit references to coordinates. Our plan is to study next the quantum dynamics dictated by the 
Hamiltonian (flOl ). on a particular surface. Such an analysis is analogous to the investigation of 
vibrational dynamics on a particular electronic potential surface, in the traditional application of 
the BO approximation. 



In steady state, the energy current operator, e.g., at the L contact, can be defined as [|44ll 

Jl = i[Hl,W], (11) 

with the expectation value 

JL{t) = Tx[JLPB{t)] = Tr[e*^^o*Jz.e-^^^°VB(0)]. (12) 

The left expression is written in the Schrodinger picture; the second is in the Heisenberg repre- 
sentation. The trace is performed over the two baths' degrees of freedom. When system-baths 
couplings, absorbed into W, are weak, the time evolution operator can be approximated by the 
first order term 

^-iHBot = e-K^L+H«)t fl-i f W{r)dA , (13) 

and the current (fT2l) reduces to 

Mt) = -i f Tr{[Ji(r), W]pLPn}dT. (14) 

Here W{t) and Jl{t) are interaction picture operators, 0{t) = e^HBtQ^-iHgt ^^^^i Hb = 
Hl + Hr. We are interested in steady-state quantities, J = Jiit ^ oo), if the limit exists. 
Expression (fT4l) can be further customized, recalling the bipartite interaction form of Vy in the 
original Hamiltonian ([T]). Then, one can formally expand W in terms of the bath operators which 
are coupled to the subsystem, B^, 

W = Y,KbBl®B], (15) 

a,b 

= E E E Kh{Bl)UB'R)ps \kp) {ms\ . 

a,b k,m p,s 

The operators Bl and Br depend on the bath coordinates, collected into Ql and Qr, respectively. 
The actual form is not important at this (formal) stage. It is specified only once particular models 
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are constructed, see e.g., Eq. (l23l) . The coefficients Aa,b absorb the subsystem parameters, the 
energies e, k and 5 in the chain model and the system-bath interaction strength X^. The powers a 
and b are positive integers. |A;) and \m) represent the many body states of the left reservoir with 
energies E^ and Em, (i.e. H^ = J2 ^k\k) {k\). Similarly, \p) and \s) are the many body states of 
the right reservoir with energies Ep and Eg. Assuming a weak system-bath coupling strength, we 
truncate W and consider only the lowest order term in BlBr, 

W ~ Ao,o + A^^^BlBr + 0{Bl) + 0{Bl). (16) 



A more general derivation in presented in Ref. uM- It can be shown that only terms containing 



products of Bl and Br add to the current, thus only the second term in Eq. (fT6l) actually matters 
for the energy current calculations. We also note that Ai i is proportional to the product Xl^r, see 
Eq. (HI). We therefore define the function T(e, k) through the relation 

Ai,i = XLXRT{e,K), (17) 



where we explicitly indicate its dependence on the subsystem parameters. Back to (1141) . employing 
Eqs. (fTTI) and (fT6l) . we obtain 

\2 /-oo 



J = '^^4^ r dt\y2xlE,r^e'^^-\BL)km{BL)mke 



(18) 



X 

p,s 

where, e.g., Z^ = J2k^~^^^'' ^^ ^^^ ^ ^^^'^ partition function; (3^ = l/T^, ks = 1, and E^m = 
Ek — Em- Time integration can be readily performed, leading to the steady state heat current 

J=^^^y udu[kL+i^)kR^iu) - kL^iu)kR+iu)]. (19) 

The excitation (+) and relaxation (— ) rate constants are given by 

Xl[{BL)km.iBL)mkfS{Ek - E^T w)^^- (20) 

I ^^ 

We have introduced here the short notation [{Bi)km{BL)mk\^ , to denote matrix elements when 
Ek > Em- Similarly, [{BL)km{BL)mkr describes the Ek < Em case. Analogous expressions hold 
for the R rates. We can also rewrite the rate constants as Fourier transforms of bath correlation 
functions 

/oo 
e^-*Tr[p,.Bi(t)i?i(0)]rft, (21) 
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satisfying detailed balance, ki+iuj) = kL-{(^)e ^^^ . Using this relation, we organize Eq. (fT9l) as 

J = V^^l!^ r ujdujkL^{uj)kR^{uj){e-^^'' - e-^«"). (22) 

This result is given in the form of a "generalized Landauer formula": The net heat current is 
given as the difference between left-moving and right-moving excitations, nevertheless, unlike 
the original Landauer formula BOll . this expression can incorporate anharmonic effects within the 
chain model, eventually absorbed into T, and nonlinear system-bath interactions, taken in by the 
rates k^±{u). We emphasize the broad status of Eq. (|22|) : It does not assume a particular structure 
for the subsystem, or a specific system-bath interaction form, 3,^, both contained inside W. It is 
valid as long as (i) there exists a timescale separation between the subsystem motion (fast) and the 
reservoirs dynamics (slow), and (ii) system-bath interactions, given in a bipartite form, are weak, 
see Eqs. ^ and (fT6l) . 

Eq. (I22I) readily reveals the dependence of the current on the subsystem parameters, thus it 
is immensely useful for exploring transport behavior. It includes a product of two terms: The 
prefactor depends on the subsystem parameters, the integral over frequencies encompasses the 
bath operators within the Fermi golden rule rates (possibly nonlinear in the bath coordinates). 
The effect of the reservoirs' temperatures is enclosed there. Since the prefactor T(e, k) is the only 
term corroborating chain parameters, by obtaining the ground state surface W [Eqs. (fT6l)-(fT7l)1. the 
scaling of the current with size and energy can be gained, without solving a dynamical problem. 

We can also regard Eq. (122)) as a generalization of the nonadiabatic transition rate, kda = 
'2iT\Vda\'^ FCW D , describing electron or energy transfer processes within a donor-bride-acceptor 
complex, to current carrying steady state situations. Here, the Franck-Condon factor FCWD ac- 
counting for the conservation of energy, depends on the temperature of the environment \JA] . In 
Eq. (I22I) this term is portrayed by the frequency integral, considering a transport process origi- 
nating from a particular state within the L bath. The second part, Vda, combines the electronic 
coupling between the donor and acceptor states. In the present work this factor is accounted for 
by the function T(e, k). 

As an example of the utility of the ETBO method to describe off-resonance conduction pro- 
cesses, consider a harmonic impurity of frequency f2, linearly coupled to two harmonic reservoirs, 

Hs = nb^b, 

H, = Y1 ""AA^^ ^- = (^' + ^)^-^- (23) 
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with B^ = ^Ablj + Kj). Here 6j, (Kj) are the creation (annihilation) operators of the mode 
j in the i' bath, b^ and b are the respective subsystem operators. Since the model is fully har- 
monic, in principle the energy current can be exactly obtained. However, this calculation requires 
some effort, and the scaling of the current with size is not easy to reveal llld. 11611 . Focusing 
on the off-resonance limit, the ETBO method can readily provide the behavior of the current 
at weak couplings. We diagonalize Hf = Hs + V and resolve the ground potential surface 
W = —^{XlBi + XrB^Y, thus extract T oc 1/fi. Relaying on the bilinear interaction form, 
the transition rates (12TI) can be obtained, A;i,+ (aj) = V y{ijj)ny{uj) , where n^^iuj) = [e^'^"^ — 1]"^ 
is the Bose-Einstein distribution function and the coefficient Tiy(uj) incorporates the system-bath 
interaction strength and the bath's density of states, assumed to be weak T^{uj) < f2. Combining 
these elements in the expression for the heat current (l22l) . we conclude that 

1 f°° 
Joe— ^/ udujTL{oj)TR{u)[nL{uj) -ur^u)]. (24) 

To be consistent with the off-resonance assumption, one should evaluate this expression at low 
temperatures T^ < fi, or impose a cutoff for the reservoirs frequencies, fi ^ coc- This result 
exposes the scaling of the current with the subsystem energy and the baths temperatures. It can be 



450. This 



shown that similar scaling holds for the spin-boson model in the off -resonance limit [|28L 
correspondence is physically correct since at low temperature an harmonic impurity behaves sim- 
ilarly to a spin impurity, as transport takes place through the lowest excitations of the subsystem. 

B. Analytic Results 

We apply the ETBO formalism on the spin-chain Hamiltonian ©-©, to obtain the energy 
current characteristics. Our objectives are (i) to resolve the behavior of the current as a function 
of chain size, and (ii) to obtain its dependence on the subsystem energetics, e and k. With this 
at hand, we can identify the dominant transport mechanism. For simplicity, we exemplify our 
analysis using the isotropic XY chain model. However, the results are applicable for other models 
including the anisotropic Heisenberg model as well, for (5 < 1, see discussion below Eq. dS]). We 
comment on this model below Eq. (l40l) . 

We recall that the basic ingredient of the ETBO formalism is the ground potential surface W, 
or its expansion, (fT6l) . Then, identifying the coefficient T(e, k), the energy and size dependent of 
the current can be captured using Eq. (I22l) . We review the elements of our model introducing a 
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more compact notation for the chain subsystem, 

Hs = eM + Kh, (25) 

with M = J2j=i ^t^J ^^'^ ^ ^^ ^^^ hopping Hamiltonian, including nearest-neighbor interactions 



ll46ll . For example, h = ^Yl {'^j'^j+i + '^J'^J+i)- ^^^ chain is connected by V^ to the thermal 



bath H^. The total Hamiltonian is given by 

H = Hf + Hl + Hji, 
Hf = Hs + V, V = Vl + Vr, (26) 

with Vl = SiBl and Vr = SnBr, Si^n oc afj^ contains subsystem operators. The energy surface 
W, the lowest eigenenergy of Hf, is accomplished through the eigenvalue equation 

Hf\go) = W\go). (27) 



Since an exact diagonalization is limited to simple models 112811 . in this work we construct W 
using time independent perturbation theory. As the unperturbed basis we utilize the subsystem 
eigenstates \n), satisfying 

Hs\n)=Er,.\n). (28) 

The system-bath interaction operator V plays the role of a perturbation. These \n) states include 
different number of excitations, demonstrated in Fig. [2l For example, for a two-qubit chain, 

Hs = e(a+ar + a^a,) + ^{a^,a^, + afa^), (29) 

we obtain the eigenfunctions and respective energies 

|0) = \U) ,^0 = 

|1) = ^(lit) - Iti)), E, = e-K 

|2) = -^(lit) + Iti)), E, = e + K 

|3) = in), E, = 2e. (30) 

The ground state is fully polarized, | H), with the two spins in their ground state. The first two 
excited states include a single excitation (a superposition, residing on the first and second sites). 
The high energy state includes two excitations. Back to the A^-site chain, the ground state energy 
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of Hf = Hs + V can be written by using time independent perturbation theory to the second order 
correction, 

W^E, + {0\V\0) + J2 ^^^J^jf . (31) 



,„ Eo — En 



The corresponding eigenfunction is 



Consider now a family of spin Hamiltonians where the ground state is fully polarized as in (|30l) . 
|0) = |i, i, •••,!)• The structure of the system-bath interaction operator, aijy, allows to connect 
this ground state only to single excitation states, the subgroup \ni) E \n), written as 

N 

\n,) = Y,C^\j). (33) 

i=i 

These are linear combinations of single excitation states, \j) = ||, |, ., tj, .., i). The eigenenergies 

of these states are En^ = e + nam with a„^ as numerical coefficients. For example, a„^ = ±1 

for the 2-qubit chain of Eq. (l30l) . Identifying the relevant states |ni), we now plug them and their 

corresponding energies into Eq. (|3TI) . The constant shift Eq was set here to zero, the second term 

vanishes. Therefore, the ground potential surface is given by 

mv\n,)f 



w ^J2 



711 



Eq — En^ 



= — 2^ . , K^ + <-'(^L) + <-'l^fi)- (^4) 

Terms which involve only B^ or B^ operators do not contribute to the current and are therefore 
ignored. Focusing on the sum, denoted by S, we simplify it recalling that n < e. We expand the 
denominator using the geometric sum formula, X]^o ^^ ~ i3^' 

"•1,9 
ni,g 

Here, the states |1) and |A^) refer to a j-type state as defined below Eq. (l33l) . containing a single 
excitation in the leftmost (1) site or in the rightmost (A^) site. The second line was derived using 
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the eigenequation for the hopping operator, h\ni) = Un^ \ni). The last line was obtained by using 
the fact that {nj\N) = and {nj\l) = for j > 2, where \nj) denotes states with j excitations 
residing on the chain. The completeness relation is also invoked, / = ^ \'n'){n\. 

We can further simplify Eq. (l35l) . We note that h is the inter-site hopping operator, and use the 
fact that it includes nearest-neighbor interactions only. This leads to {1\ h'' \N) = if g < (A^ — 1). 
Therefore, the leading term of the q expansion must include A^ — 1 operators for transferring an 
excitation from the first unit of the chain to the last one, 

1 / K\^-l 



W^-- (^--j BlBr [{l\ h""-' \N) + {N\ /i^-i |1) 



(36) 



The square brackets yield a numerical factor. We conclude that the ground state potential follows 
a simple form 

W^r{e,^)BLBn, (37) 

with 

r(e,«:) = --(--) . (38) 

We now go back to the energy current (l22l) . denoting by F(Ti, Tr) the contribution that depends 
on the reservoirs' temperatures, 

J = ^[-) F{Tl,Tr). (39) 

One can also express the prefactor by a decaying exponent, 

r(e, k) oc e""^, a = -21n(fi:/e). (40) 

Eq. (l39l) describes an exponential decay of the energy current with distance, a coherent- 



superexchange dynamics Ill5|,l35l]. The physical picture exposed is that low frequency (reservoirs) 
modes are being coherently exchanged, without the actual excitation of the (off-resonance) modes 
of the chain. The intermediating structure-chain therefore serves as a mediating medium, allow- 
ing for through-bond couplings. This expression is the analog of the McConnell super-exchange 
result, describing deep electron tunneling in tight binding models [|35l]. 

The derivation above is applicable for models more general than the Hamiltonian ([I])-©. For 
example, we may modify the spin chain and add an interaction term 5(7|cr|^]^; The exponential 
result still holds as we show next. Overall, the analysis relays on the following ingredients: (i) 
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FIG. 4: Potential energy surfaces for iV = 2 (top) and iV = 4 (bottom), e = 3, k = 0.2, 5 = 0, Br = 0.1. 




FIG. 5: Contour plot of the four potential energy surfaces for N = 2. Panels (a) to (d) present the potentials 
from the ground state upwai^d, e = 3, k = 0.2, and 5 = 0. 

the ground state of Hs is a fully polarized state, i.e., there are no excitations on the bridge in the 
absence of the interaction with the reservoirs, (ii) system-bath interactions involve the generation 
of a single excitation on the chain boundary sites, and (iii) the chain units are weakly connected 
through nearest-neighbor couplings, small relative to the onsite gap, k <^ e. Under these condi- 
tions, combined with the perquisites for the validity of the ETBO approach (off-resonance condi- 
tion and weak system-bath couplings), transport dynamics reflects the superexchange mechanism, 
Eq. (EOl). 

The above analysis holds, under some conditions, for describing the dynamics of the anisotropic 
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Heisenberg chain in the off-resonance regime. For a 2-qubit model 

Hs = € [aia^ + a^a^) + ^ {a^a^ + alal + 5ala'^) . (41) 



We repeat the derivation, Eqs. (I34l)-(l35l). and resolve the ground potential 

{Bl-BrY , {BL + Bnf 



W 



^ [(e - k6) Bl - 2kBlBr + (e - k6) BI] 



k5Y — K 



2 



H _K^ K^ 






e2 e3 ' e^ 



Sz.5ii + 0(i?i) + 0(i?i). (42) 



The last line was derived under the weak exchange assumption, e :» k. This result agrees with 
the behavior predicted in Eqs. (I37l)-(l38l) when 5 = 0. We also note that the exchange anisotropy 
parameter 5 affects W to higher order in n/e. Thus, to the lowest order in n/e the energy current 

2 

satisfies J{N = 2) oc ^F{TL,Tf{), irrespective of the details of the spin model. This behavior 
prevails for longer chains as well: The onset of S provides only higher order corrections to the 
off -resonant energy current (l39l) when 6 is small. More precisely, the results hold as long as the 
spectrum maintains its distinct band structure, see Fig. [3l the (single) excitation energies could be 
still approximated by -E^^ ~ e + fi:a„^ , and the ground state is fully polarized. 

We now comment on the usefulness of the ETBO method to describe transient effects. While 



the approach has been formulated for treating non-equilibrium steady state situations [12811 . one 
could also rewrite it to describe the transient dynamics of a donor excitation, transferred to an 
acceptor sidegroup through a bridging backbone. In the context of electron transfer, twp distinct 
quantities, the electrical conduction and the electron transfer rate, were shown to be linearl y re - 



lated [|33l]. Similar correspondence should arise in the context of excitation energy transfer [|15|1 . 
comparing the steady state energy current at very small temperature bias and the excitation transfer 
rate. We thus argue that the thermal conductance, obtained as limAr=o J/^T, isproportional to 



the excitation energy transfer rate detected in donor-bridge- acceptor complexes 113 311 . 



C. Numerical Simulations 

We support our analysis by an exact numerical diagonalization of Hf, to obtain the set of po- 
tential surfaces Wn- We recall that the potential surfaces Wn{Q), with Q enclosing the bath coor- 
dinates coupled to the system, are the analogs of the electronic surfaces in the context of molecular 
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FIG. 6: The function V = BLBiifs{e, k) for N = 1,2, 3, 4. The bath coordinates are fixed at Bl = Bji 
0.05. The exchange parameters are k = 0.1, 6 = 0. 




FIG. 7: The function V = BiB^fsie, k) for A^ = 2, 3, 4. The bath coordinates are fixed at Bl = Br = 
0.05. The exchange parameters ai^e e = 4, 5 = 0. 

structure, generated by varying the (slow) nuclear coordinates. Here, in a similar fashion we treat 
the bath coordinates Q as parameters, overall treating Bl{Q) and Br{Q) as parameters. The 
model consists the isotropic XY spin chain coupled at the boundaries to two thermal reservoirs, 
Eqs. (HI)-©. First, in support of the adiabatic approximation we show in Fig. |4]that the ground 
potential energy surface W is separated by a substantial gap from other states. The x-axis is the 
bath coordinate B. Practically, the data was generated by fixing 5/j and parametrically modifying 
the coordinate Bl. The ground potential surface W lies around —Ne/2. It is separated by ~ e 
from the higher excitation states. This observation consistently supports the BO scheme. While in 
Fig. |4]the potentials seem flat due to the scale used, in Fig. [5]we explicitly present a contour plot 
of the four potential surfaces Wn for A^ = 2, to demonstrate their variance with the thermal bath 
coordinate By. We now select W, the ground potential, and postulate that it fulfills Eq. (fT6l ). or 
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more generally, 

W = fo{e) + h{Ble,^) + fn{Ble,K) 

+ BLBRfs{e,^) + 0{BlBl), (43) 

with the unknown functions /. For resolving the part in W responsible for transport, we define an 
auxiliary function 

W' = W-W{k = 0). (44) 

The difference V = W'- W'{Bl = 0) - W'{Br = 0) should isolate the (fourth) term in Eq. dH, 
the term which contributes to the energy current in the lowest order of the system-bath interaction 
strength. Figs. [6l|7] display this function for chains with A^ =1,2,3, and 4 units We conclude that 
the exponential law, Eq. (l38l) . is indeed satisfied. Fig. [6] verifies the exponential dependence on e, 
whereas Fig. |7]proves the same for the intersite coupling k. 

IV. RESONANCE TRANSPORT: MASTER EQUATION FORMALISM 

A. Method 

Our objective here is to simulate resonant energy transfer across isotropic XY spin chains, 
assuming that the bath populated modes overlap with the subsystem gaps. The dynamics is inves- 
tigated using the Born-Markov approximation, a second order perturbation theory scheme which 



invokes the Markov approximation ll47ll . Furthermore, using the secular approximation (SA), the 



diagonal and nondiagonal elements of the density matrix are separated. This scheme results in a 



markovian quantum master equation. The method is detailed in Ref. 114311 where it was applied 
onto an impurity single-site model. Here, we generalize this treatment for studying the energy 
current behavior in an extended system. Comments about the validity of this approach, to describe 
energy transfer processes in spin chains, are included below Eq. (|52l) . 
We begin by diagonalizing the subsystem Hamiltonian 

Hs = LHsLK (45) 

L is a unitary matrix which diagonalizes Hs- As before, we denote the resulting eigenspectrum 
by \n) with En- The subsystem operators coupled to the bath are transformed to the new basis. 

Si = L^a^L, Sn = L^a%L. (46) 
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The operators 5* can be formally presented by their matrix elements as 

Si = y^5'i,mn|m)(n| 
nm 

Sn = '^SN,mn\m){n\. (47) 

nm 

Since the system is uniform it can be shown that l^i^mnp = ISj^^rnnl"^ = |'S',„„p, where the site 
index is neglected. The total Hamiltonian in the subsystem basis is given by 

H = Hs + \l~SiBl + Xb~SnBr + Hl + Hr. (48) 



Under the Born-Markov scheme [|47h . accompanied by the SA, the probability to occupy the n 



subsystem state can be described by a first order differential equation 

-< n = / ^ |*-'mn| -T m (, * j ft;„^_^„ — rn\t) / ^ \Jm.n\ l^n^rn- V^") 



v,m 



The transition rate constants satisfy [|47ll 



/oo 
dte-'^-'i:iB[B,{t)BMl (50) 

-oo 

where Enm = En — Em and the operators are written in the interaction representation, B^(t) = 
^iH^tj^^^-iH^t^ rpj^g trace is performed over the L and R bath states. In steady state, the set (|49l) 
reduces into a linear set of equations. Complemented by the conservation of the total probability, 
^ P„ = 1, we can numerically obtain the steady state occupation probabilities at each state. The 
energy current, at the level of the Bom-Markov approximation, is given by [|43|] 



n,m 

It can be readily calculated once the steady state population and rate constants are known. At this 
stage one should choose a particular form for the bath operators coupled to the subsystem. For 
example, selecting the displacement operators ll43n . the rate constants reduce to (m > n) 



C^„ = T{Emn)[nuiEmn) + 1], (52) 

with T{u) = 2TrXl ^ . 5{u — Uj). In practice, we take F as a constant, independent of frequency, 
identical at the two contacts. The function n^(u;) = [e'^/-^" — 1]"^ is the Bose-Einstein occupation 
factor. 
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The authors of Ref. 14811 questioned the validity of a related approach, the Redfield equation, 
derived in the chain-local basis, for describing the dynamics of several spin chain systems, hi 
particular, under the secular approximation, zero energy current was obtained in nonequilibrium 



situations ll48h . Inconsistencies of the Redfield equation, unable to properly reproduce equilibrium 



and nonequilibrium dynamics, were noted in the past in the context of electron transfer processes. 



see e.g., Ref. ll49ll . There, it was argued that working in the subsystem eigenbasis should lead 



to a proper equilibration process and to the correct nonequilibrium dynamics. In view of the 



zero-current at finite bias anomaly ll48n . we work here in the chain eigenbasis, indeed naturally 



eliminating such a nonphysical behavior 



The authors of Ref. [14811 further traced the nonphysical dynamics within the Redfield approach 
to the inconsistency of the secular approximation, when applied onto the chain model. It is argued, 
that this approximation, resulting in the separation between the diagonal and nondiagonal terms 
of the reduced density matrix, relays on the assumption that differences between the subsystem 
energies are large compared to the subsystem relaxation rate constants. However, in the chain 
model differences between energy states within each band diminish for long chains, thus one 
should carefully review the SA, as we do next. 

The eigenspectrum of the isotropic XY spin chain was presented in Fig. [2l We recall that for 
e ^ K the subsystem's energies are grouped into manifolds, each containing a particular number 
of excitations. It should be noted that the gaps between bands are preserved, order of e, even 
for long chains. We argue next that even though within each manifold the states become quite 
dense, one could obtain the correct dynamics of the isotropic XY model using standard quantum 
master equation approaches, stating the SA, as long as the gaps between different bands are main- 
tained. The reasoning is that once we work in the chain-diagonal basis, the equation of motion 
for the density matrix (before the SA) connects only states which differ by exactly one excitation 
through bath excitation and relaxation processes. Rephrased, states within the same manifold are 
not directly linked, only through higher-order bath correlation functions. Thus, within the Bom 
approximation, energy differences that come into play within the density matrix equations are al- 
ways order of the gap e. Since we pick small relaxation rate constants F < e, we conclude that the 
SA is consistent in the present setup. This argument does not hold for the Heisenberg model, as 
the excitation gaps rapidly disappear with increasing size, see Fig. [3l 
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FIG. 8: Left panel: Energy current as a function of size for the isotropic XY chain in the resonant Umit, 
e = 0.5, 1, and 2, bottom to top. Right panel: Energy current as a function of spin gap for A^ = 6. Other 
parameters are k = 0.1, T = 0.01, T^ = 4 and Tr = 2. 
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FIG. 9: Energy current as a function of bridge size assuming a superexchange (full), or a ballistic (dashed) 
mechanism. The dotted line is the total current. Other parameters are e = 2, k = 0.2, F = 0.05, Tl = 0.2 
andTR = 0.1. 

B. Numerical Simulations 

The steady state dynamics of the isotropic XY model in the resonant regime is presented in Fig. 
[8l The left panel displays the energy current as a function of chain size. We note that the current 
scales as J oc N^, indicating on a ballistic energy transfer mechanism [|48|] . The right panel of Fig. 
[8] presents the behavior of the current as a function of spin gap e. At high temperatures, T^ > e, 
the current follows J oc e, as expected for a ballistic motion. For large e, beyond the reservoirs 
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temperatures, the current declines since many (high energy) subsystem modes cannot participate 
in the transport process any longer as bath modes matching the subsystem frequencies are not 
significantly populated. Qualitatively, one may suggest that J oc e/(e^ + a^), where a is large at 
high temperatures. 

As discussed above, the master equation method followed here cannot be utilized once the large 
gaps in the band structure close, see Fig. ^ N = 10. Therefore, we cannot faithfully describe 
here the role of the anisotropy exchange parameter 5 on the dynamics. Using a Redfield type 
approach without invoking the SA, it can be shown that for larg e enough 5, instead of the (resonant) 
ballistic dynamics, heat propagates in a diffusive manner [21]. Therefore, while the off -resonance 
superexchange dynamics, relaying on the bridge as a mediating medium, does not depend on the 
fine details of the chain Hamiltonian, in the resonance regime transport characteristics crucially 
depend on the details of the chain structure. 

V, CONCLUSIONS 

We studied the energy transfer behavior in homogeneous linear spin chain models coupled 
at the two ends to thermal reservoirs in two opposite limit: in the off-resonance and resonance 
cases. In the off-resonance limit the dynamics was investigated by adopting the recently developed 



ETBO method [|28|] . The combination of analytic manipulations and numerical simulations con- 
firmed that the energy current exponentially decreased with distance, an indication of a coherent- 
superexchange transport mechanism. This behavior is generic, irrespective of the details of the 
chain model. In the resonant regime a standard master equation method was used, specifically 
demonstrating that the energy dynamics in the isotropic XY chain model is ballistic, as the current 
does not depend on the system size. 

We separately presented theories for describing off -resonance and resonance energy transmis- 
sion, with the bridge modes located either above or in resonance with the reservoirs populated 
modes. A complete theory for describing, on the same footing, these two limits could be based on 
a surface hopping approach lISOll . or relaying on a nonmarkovian master equations for describing 
the chain dynamics ll47n . Here we demonstrate the crossover between the superexchange behav- 
ior and the resonant dynamics by showing, on the same plot, the deep-tunneling energy current, 
the ballistic component, and the total current, as a function of bridge length, see Fig. [9l Data 
was generated for the isotropic XY chain connected to ohmic-bosonic reservoirs maintained at 
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low temperatures. The superexchange behavior was simulated by adopting Eq. (1221) . The ballis- 
tic component was gained using the method explained in Sec. IV. We find that for short chains 
the coherent- superexchange contribution, resulting from the transmission of low frequency modes 
across the bridge, dominates the current. In contrast, for long chains resonant conduction is more 
significant, though the population of bath modes matching the system gaps is small at low temper- 
atures. The turnover between the tunneling dynamics and the resonant behavior occurs between 
N = 1 to 2 for a broad range of parameters, e = 1 - 2, k = 0.05 - 0.2, T = 0.01 - 0.05, 
T,^ ~ 0.1 — 0.5 (dimensionless units of energy, h = 1). This observation lies in general agreement 
with recent experiments of triplet energy transfer on vr-stacked molecules, demonstrating that the 



turnover between tunneling and (resonant) diffusive mechanisms occurs between A^ =1 to 2 [|40l1 . 
We expect that the Heisenberg model will similarly show a turnover between the superexchange 
mechanism and the diffusive (hopping) dynamics around similar bridge sizes. 

While the present analysis was mainly carried out adopting the isotropic XY chain as the bridg- 
ing object, the results of the ETBO method hold for the anisotropic Heisenberg chain and other 
similar variants, as long as gaps between different excitation manifolds are larger than energy 
differences within each band. Furthermore, the total Hamiltonian, combining the reservoirs and 
(nonlinear) system-bath interactions, cannot be generally mapped onto a noninteracting fermion 



model 114 in . 



The energy tunneling- superexchange behavior observed in the off-resonance regime has been 



discussed before in the context of excitation energy transfer [|15|] . Here it is rigorously obtained 
in a first principle derivation, relaying on the timescale separation between subsystem dynamics 
and the baths' motion, irrespective of the details on the chain spectrum, the reservoir realization, 
and system-baths interaction form. We expect this general behavior to show itself in numerous 
systems, including organic and biological structures, exploring electronic B91 14011 and vibrational 
M\ energy transmission. 
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